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We consider the problem of digital halftoning from the view point of statistical mechanics. The 
digital halftoning is a sort of image processing, namely, representing each grayscale in terms of black 
and white binary dots. The digital halftoning is achieved by making use of the threshold mask, 
namely, for each pixel, the halftoned binary pixel is determined as black if the original grayscale 
pixel is greater than or equal to the mask value and is determined as white vice versa. To determine 
the optimal value of the mask on each pixel for a given original grayscale image, we first assume that 
the human-eyes might recognize the black and white binary halftoned image as the corresponding 
grayscale one by linear filters. The Hamiltonian is constructed as a distance between the original 
and the recognized images which is written in terms of the threshold mask. We are confirmed that 
the system described by the Hamiltonian is regarded as a kind of antiferromagnetic Ising model 
with quenched disorders. By searching the ground state of the Hamiltonian, we obtain the optimal 
threshold mask and the resultant halftoned binary dots simultaneously. From the power-spectrum 
analysis, we find that the resultant binary dots image is physiologically plausible from the view 
point of human-eyes modulation properties. We also propose a theoretical framework to investigate 
statistical performance of inverse digital halftoning, that is, the inverse process of halftoning. The 
inverse-halftoning is regarded as a special example of image restoration in which one should infer the 
original grayscale image from the less informative black and white binary dots. From the Bayesian 
inference view point, we rigorously show that the Bayes-optimal inverse- halftoning is achieved on a 
specific condition which is very similar to the so-called Nishimori line in the research field of spin 
glasses. Finally we show that both halftoning and the inverse-halftoning processes are unified under 
a single Hamiltonian, namely, it is possible for us to obtain the threshold mask, the halftoned and 
inverse-halftoned images simultaneously by finding the ground state of the spin systems. 
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I. INTRODUCTION 

Recently, a lot of problems of information science and technology have been investigated by several useful tools 
developed in the research field of statistical mechanics of spin glasses [3-[3]- Statistical mechanics of information is 
now widely spreading in various subjects such as neural networks Q, error-correcting codes 043; CDMA multi-user 
demodulator [8j, image processing, etc HQ. Especially, image restoration by making use of a graphical model of 
Markov random fields has been investigated extensively from both analytical and numerical point of views @. In 
the research field of image processing, digital halftoning instead of image restoration, which is denned as a process 
of generating a pattern of pixels with limited number of colors, especially converting a grayscale image into the 
binary black and white picture, has been widely used in various practical situations in media, such as the printing of 
newspapers, fax machines and so forth [Tol fTlj. 

To achieve the digital halftoning, one needs the strategy to arrange the geometrical-combination of black and white 
pixels so as to make human-eyes to have a kind of optical illusion. Namely, the halftoning relies on the fact that 
the human-eyes act as a spatial low-pass filter and can not recognize any specified structure in the part of image 
dominated by high frequency components. 
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Actually, to justify this fact, several authors proposed or estimated the so-called the contrast sensitivity function, 
or the modulation transfer function of human visual systems independently. For instance, Analoui and Allebach (l2l | 
introduced the contrast sensitivity function having the following form: 

F(f p ) = fc{exp(-27ra/ p )-exp(-27r/3/ p )} (1) 

where f p denotes frequency and a and (3 are constants setting to 0.012 and 0.046 respectively. They set the scale 
parameter k to satisfy the condition max /p F(f p ) = 1, namely, k = [{aJB)- a/{ - a - l3) - {a/ '{})-&/ {"-P)]- 1 ~ 2.173938. 
As another model of the human visual system, Mannos and Sakrison [13f estimated the function as 

F{f p ) = 2.6(0.0192 + 0.114/,) exp [-(0.114/p) 11 ] . (2) 

We plot the shape of these contrast sensitivity functions ([1} and ^ in FIGfT] From this figure, we find that in 




FIG. 1: Human-eye modulation transfer functions estimated by Analoui and Allebach 1 1 21 (the solid line), and Mannos and Sakrison 
[l3l](the broken line). 

high-frequency regime / p > 1, the both contrast sensitivity functions decrease to zero, which means that human eyes 
cannot recognize high-frequency components in images. In other words, the halftone algorithm should be constructed 
such that the halftoned binary dots contain relatively high-frequency components to describe the original grayscale 
revels. This fact is an important guide for us to consider halftoning algorithms from the view point of physiology. 

Up to now, to achieve fine qualities of the digital halftoning, a lot of techniques, for instance, clustered- dot ordered 
dtiher method, threshold mask method 0, EH and blue noise mask method (l4l Ila]. error- diffusion method [l6j . etc. 
have been proposed and developed by many researchers. However, no attempt has been done to formulate the 
important problem from the view point of statistical mechanics of information. It is not difficult for us to assume 
that each pixel for the halftone and the grayscale images is represented by Ising spin and Q-Ising spin (or Potts spin) , 
respectively. Therefore, it seems that statistical mechanical approach is very useful for the problem. 

On the other hand, to retrieve the original grayscale image from the halftone binary dots, the so-called inverse- 
halftoning has been also developed especially in scanner technology [13, [3 ■ As such an inverse halftoning strategy, 
the conventional smoothing filters have been used widely. However, it is worth while for us to look for the alternative 
or reconsider the problem from the view point of statistical mechanics of disordered spin systems. 

In this paper, we propose a statistical-mechanical modeling of both digital halftoning and the inverse process, 
namely, inverse-halftoning. To attempt to generate the halftoned binary dots, we first introduce a threshold mask 
matrix with the same size as the original image. For the threshold mask, each pixel of the halftoned image is 
determined as a black pixel if the value of the original grayscale pixel is greater than the component of the threshold 
mask matrix and is determined as a white pixel vice versa. Then, we assume that human eyes might recognize the 
original grayscale image by a linear filter. By taking into account that the distance between the original and the 
recognized images is written in terms of the threshold mask, we naturally introduce a Hamiltonian to be minimized 
to determine the optimal threshold mask as the ground state. We explicitly show that the system described by the 
Hamiltonian is a variant of the antiferromagnetic Ising model with disorders. For a demonstration of our method 
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to generate the optimal threshold mask, computer simulations of simulated annealing [19|, l20( is carried out. To 
investigate the statistics of the threshold mask image and the halftoned image, we evaluate the power-spectrum 
numerically. 

On the other hand, in order to investigate statistical performance of the inverse digital halftoning, in our previous 
study Jj) , we introduced the standard regularization theory for a Markov random-fields model to represent original 
grayscale image and construct the inverse process as a kind of image reconstruction of the grayscale image from a given 
halftone binary dots. Then, we evaluated the statistical performance by making use of Markov chain Monte Carlo 
(MCMC) simulations and analysis of the infinite-range mean-field model. We also investigated the inverse-halftoning 
process as a dynamics of disordered spin systems. In this paper, we reconsider the inverse halftoning problem from 
the Bayesian inference point of view. We show that the Bayes-optimal inverse-halftoning that minimizes the mean 
square error is achieved on a specific condition which is similar to the so-called Nishimori line [2lj in the research 
field of spin glasses. Finally we discuss the possibility that both halftoning and the inverse-halftoning processes are 
unified under a single Hamiltonian, namely, we argue whether it is possible for us to obtain the threshold mask, the 
halftoned and inverse-halftoned images simultaneously by finding the ground state of the spin systems or not. 

This paper is organized as follows. In the next Sec. UTJ we explain our procedure to generate the optimal threshold 
mask for a given original grayscale image. In this section, the Hamiltonian is introduced as a function of the threshold 
mask under the assumption that the human-eyes might recognize the halftoned image as the corresponding original 
grayscale version by a linear filter. In Sec. IIII1 we solve the combinatorial optimization problem described by the 
Hamiltonian by making use of the simulated annealing. We show the resulting halftoned image for a standard grayscale 
image. We find that our algorithm generates the binary dots which induce the illusion of a continuous-tone image. 
We also investigate the power-spectrum statistics of the threshold mask array and the halftoned image. In Sec lIVl 
we formulate our inverse-halftoning procedure based on the MPM estimation of Bayesian statistics. We show clearly 
that the Bayes-optimal inverse-halftoning that minimizes the mean-square error is achieved on a specific condition. 
In Sec. El we discuss the possibility that both halftoning and the inverse-halftoning processes are unified under a 
single Hamiltonian. The last section is summary. 

II. STATISTICAL-MECHANICAL MODELING OF HALFTONING 

We first deal with the 'forward problem' of the digital halftoning. In this section, we first explain the procedure to 
generate the binary dots to represent the original grayscale image via what we call threshold dither method [lol [Tl| . 
Then, we explain how one formulates the halftoning processes as a problem of disordered spin systems and why 
statistical mechanics of information is useful for the problem. 

A. Halftoning by threshold dither method 

Let us first define the original grayscale image located on the square lattice with size L\ x L 2 by 

g = {g x , y G 0, 1, • • • , Q - l\x = 1, 2, • • • , L u y = 1, 2, • • • , L 2 } . (3) 

To generate the binary dots to represent the grayscale image g, we introduce the threshold mask matrix with the size 
h x l 2 {h < L\M < L 2 ): 

t = {^eO,V - ,Q-l\x = l,2,--- ,Zi,i/ = l,2,... ,h] (4) 

We should keep in mind that each component of the matrix t x ,y should be satisfied the following periodic boundary 
condition: 

tx, y =t x ±h, y = t x>y ±i a , L\ = (mod h), L 2 = (mod l 2 ). (5) 
Then, the halftoned image is defined by 

h = {h x , y €0,l\x= 1,2,- ■■ ,L u y = 1,2,- ■■ ,L 2 } (6) 
and each pixel in the h is calculated as 

^x^y hx,y ®(.9x,y ^x.y)i (7) 
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where the step function Q(x) is defined conventionally as 

Q(x) = 



1 [x > 0) 
(x < 0) 



(8) 



For this setup of the halftoning procedure, a pixel of the resulting halftoned binary image is set to one if the pixel of 
the original grayscale image is greater than or equal to the corresponding pixel of the threshold mask; otherwise the 
pixel is set to zero. 

For instance, the digital halftoning is achieved by threshold mask of Bayers' type for Q = 16 with size iixl 2 = 4x4: 
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(9) 



Namely, for a part with size 4 x 4 of the original image g, say, for 
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we have the corresponding block of the halftone image h as 
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In FIG. [U we show a typical example of the original grayscale image (256 grayscale levels: left panel), the image 
with reduced grayscales Q — 16: center panel) and the resulting halftone image obtained by the Bayer's threshold 
mask (|9|). The size of the image is L\ x L2 = 256 x 256 and the number of the grayscale levels is 256. Here we first 
converted the original 256-grayscales g image to the reduced |N(x, y)\ — Q(< 256)-grayscale image g by the following 
transform: 



C{ 
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x : Q) = k \ 9 



■ r ~{k 
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256. 
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(12) 



and used V X)JI g x , y — C{g x y : Q) as the original pixel in ([7]). From the right panel of this figure, we find that the 
grayscale-levels of the image in the middle panel is described by black and white binary dots. To see the detail of 
the binary dots which represent the grayscale levels in the original image, in FIG. |31 we show the result focusing on 
a tiny part (the 'right eye' of the woman) of the same image in FIG. [2] 

Therefore, once we obtain the appropriately threshold matrix t, the halftone binary dots h is determined uniquely, 
and for this reason, the quality of the halftoning is dependent on how we choose the threshold mask appropriately. In 
this paper, we consider the case of the threshold mask having the same size as that of the original grayscale image, that 
is, the case of l\ = L\, I2 = L2, and propose an algorithm to obtain the optimal threshold mask as lowest energy states 
of the Hamiltonian which is a square distance between the original and the recognized images. We shall introduce it 
in the next subsection. 



B. A recognition model of human vision 

As we explained in the previous subsection, our main problem is now to determine the threshold mask t efficiently. 
For this purpose, we might assume that human eyes might retrieve the grayscale image, let us call, the recognized 
image: 

s = {s x . y G 0, 1, • • • , Q - l\x = 1, 2, • • • , L lt y = 1, 2, • • • , L 2 ] (13) 
from the halftoned image h by making use of the following linear filter: 

V ' x,y ^x,y — ^ ^Vx~i,y—j (44) 
i,jeN(x,y) 
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FIG. 2: An example of digital halftoning by using the Bayers' matrix ^S). The 256-grayscale image (left) and the original images in which 
the grayscale levels are reduced from 256 to 16 (center). The right panel shows the halftoned image obtained by the dither method of (5}. 
We set Q = 16 and L\ x L 2 = 256 x 256. 





FIG. 3: The result focusing on a tiny part (the 'right eye' of the woman) of the same image as in FIG. [2] 



where N(x, y) denotes the nearest neighboring pixels of the pixel located at (x, y) and the point (x, y) itself. Apparently, 
for a two-dimensional square lattice, the ingredients of N(x, y) are (x, y), (x — 1, y + 1), (x, y + 1), (x + 1, y + 1), (x — 
(x + 1, y), (x — l,y — 1), (x,y — 1), (x + 1, y — 1) and the size is |N(x, y)\ — 9. The above choice for the recognition 
model comes from our assumption that human-eyes might recognize a part of halftone pictures more dark when the 
density of black pixels located on the part is relatively high. The justification of this simple assumption could be 
checked by comparing the recognized image s, which is calculated in terms of the resultant halftone binary dots via 
(TH)) , with the original grayscale image g. We shall discuss the results in Sec. IHII 

Obviously, for the simplest choice, we might set the weight W x>y appearing in (fT4|) as W x>y = 1. We should 
notice that for this choice, the spatial structure of the black dots in N(x, y) is not taken into account and the number 
|N(x,y)| itself determines the grayscale level of the corresponding pixel s x . y in the recognized image s. Therefore, 
the choice of the effective grayscale level can be controlled by choosing the number of the nearest neighboring pixels 
|N(x,y)|. Taking into account this limitation, we first reduce 256-grayscale levels to Q{< 256) by (fT2|) and regard the 
Q-grayscale image as the g. 

We also might utilize the other choice such as the following two-dimensional Gaussian- type: 
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(15) 



For this choice, not only the number |N(x, j/)| but also the spatial structure determines the grayscale levels in the 
recognized image. However, we naturally assume that the spatial structure does not affect the quality of halftoning 
if the size of the window |N(x, y)\ is small enough in comparison with the image size L 2 . 

In this paper, we use the following specific choice of the weight W x -i,x-j with size |N(x, y)\ — 5 x 5 = 25 as 



(W x - 2 .y+2 

W x - 2 ,y+l W x ^!^y +1 W X: y +1 W x+ 

W X - 2 ,y W X -l,y W X: y W X+1 ,y 

W X -2, V -1 W x -l,y-l W x ,y-1 W x+ l >y 

V W x -2,y-2 W x -l, y -2 W x , y -2 W X+ l,y- 2 



w x W x+1 ,y+2 W x +2,y+2 \ 

Wx,v+1 "x+l,y+l W x +2,y+l 

" ■ W x+2 ,y 



i w 3 



x+2,y-l 



(I 



x+2,y 



-2) 



5 15 5 5 \ 

10 15 10 5 ' 

15 15 35 15 15 

5 10 15 10 5 

5 5 15 5 5 



(16) 



6 



We should keep in mind that each recognized pixel s XtV takes the minimum s Xi y = and the maximum s XiV — 255 from 
the definition ([Til) . Hence, each pixel in the recognized image s calculated by (ITT)) represents 'defective' 256-grayscale 
levels, that is, V XtV s x , y e {0, 5, f 0, • • • , 255}. 

Then, we might choose a strategy to determine the threshold mask t that minimizes the square distance between 
the original and the recognized images. Namely, we minimize the energy \\g — s\\ 2 = J2 x Li (9x,y ~ s x, y ) 2 - We 

should notice that from the relation (fTT)) . the above distance ||g — s|| 2 is written in terms of the threshold mask t for 
a given original image g. By substituting (fTT)) into the distance \\g — s|| 2 , we obtain the Hamiltonian of the system 
as a function of the threshold mask t. Thus, our digital halftoning is now reduced to a combinatorial optimization 
problem of the following Hamiltonian: 

Li L 2 ( ) 2 

= E E \ - E w»-i,v-iQ(fl^ - \ ( 17 ) 

x=ly=l [ i,j£N(x,y) J 

The above Hamiltonian H(t\g) is a starting point of statistical-mechanical modeling of the digital halftoning. We 
should bear in mind that in the above Hamiltonian, t are dynamical variables, whereas, g are quenched disorders. 
As we have the Hamiltonian of the system, one can generate statistical ensembles according to standard statistical 
mechanics. In its thermal equilibrium at temperature each possible microscopic state t obeys the following 

Boltzmann-Gibbs distribution: 

m\9) = l L eM oT l9)] (18) 

One of the simplest uses of the Hamiltonian is to regard its minimum energy state as the optimal threshold mask 
t = {t X; y € 0, l\x = 1, 2, • • • , L\, y — 1, 2, • • • , L 2 }. Namely, we have 

i = aigm&x Pp(t\g), (19) 
which is rewritten in terms of the statistical-mechanical terminology as 

Li L 2 Q-l 

y x , y L, y = Hm H II E ^M*\g)- ( 2 °) 

Calculating the Q ilxL2 -sums: Jlsii 11^=1 Y^t'^oi' " ) an( l t & king the zero-temperature limit are performed by 
Gibbs sampler with temperature annealing procedure during the Monte Carlo steps. The solution might be a good 
candidate for the suitable threshold mask. 

We should mention that our procedure is categorized in the method 'model-based halftoning' which means that 
the threshold mask is dependent on the original image. In engineering, there exist several such algorithms [22|-|2J]. 
However, we should stress that these studies are completely different from our statistical-mechanical modeling. 



C. The corresponding spin system 

It is helpful for us to consider the meaning of the Hamiltonian (| 1 7[) from the view point of disordered spin systems. 
To see the relationship between halftoning processes and the corresponding disordered spin system, we expand the 
square of the Hamiltonian (|17j) . Then, we have 

n{t\g) = - EE E E j nM s i,3 s k,i + EE 5a ^ E h ^ s ^ ( 21 ) 

x=ly=l (i,j)/(fc,i)eN(a:,y) (fc,i)£N(x,y) x =ly=l i,j£N(x,y) 

where we defined hij = W x -i, y -j > 0, J ih u = W x -j, y -jW x -k,y-i > and Sij = sgn(g y - ty) e {-1, +1}. It should 
be noted that we canceled a constant term Y^, x Li T^yLi{9x.y) 2 and used the relation Q(x) = (1 + sgn(x))/2. 

Thus, our system is nothing but an Ising model defined on the two-dimensional square lattice with antiferromagnetic 
interactions with random field g on pixel. However, the ground state which minimizes the above Hamiltonian with 
respect to t is complicated due to the quenched variables as a given original image g appearing in the argument of 
sign function such as sgn^y — ty). The lowest energy state in the dynamical variable space of t might be highly 
degenerated and it means that we can use various techniques developed in the field of statistical physics of disordered 
spin systems to obtain the lowest energy state and to investigate the lowest energy properties. In the next section, we 
show some demonstrations to make the halftoned binary dots that imitate the original grayscale image by minimizing 
the Hamiltonian (|17[) via simulated annealing. 
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III. NUMERICAL EXPERIMENTS 

In this section, we show some demonstrations to generate the optimal threshold mask and the resulting halftoned 
image for a given standard grayscale image. We shall demonstrate our algorithm for one of the well-known standard 
images shown in FIGH] (upper left). For the original image g with the size L\ x L 2 = 400 x 400 as input data, we 
construct the Hamiltonian H(t\g) and carry out the calculation of the threshold mask t given by (|20)l via simulated 
annealing with the temperature schedule j3 ~ \fi. For the solution of the threshold mask t, we make the halftoned 
image h according to (j6|). 

In FIG. 21 we show the results. From the structure of the binary dots (see FIG. H] (upper right)), we find that the 




FIG. 4: From the upper left to the lower right, the original image of 256-grayscale levels, the halftoned image, the threshold mask, and 
the recognized image obtained by our algorithm. The size of the images is L\ X L2 = 400 X 400. 

resulting halftoned image looks well globally, however, locally it contains some curious clusters of vortex configurations 
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(see FIG. [5]). The similar phenomena are generally observed in the halftoned binary dots image via the error- diffusion 
method. 




FIG. 5: A part (a local structure) of the resulting halftone image (the 'mouth' of the woman shown in FIG. [4] (upper right)). Curious 
clusters of vortex configurations are observed. 



In the lower right panel of FIG. SI we show the recognized image s. We find that the recognized image is quite 
similar to the original grayscale image and this result is a justification for our assumption to construct the recognition 
model (|14|) . From these pictures, we might conclude that our algorithm works very well and the resulting halftoned 
binary dots look fine to represent the original grayscale levels. However, it is quite important for us to evaluate the 
performance quantitatively. In following , we evaluate the performance of the halftoning quantitatively. 

In FIG. [6] (right), we show the histograms of grayscale levels for the original and recognition images. We find that 
the both shapes of the histograms have a similarity, however, there is a gap between the histograms in their height. 
Apparently, if the perfect minimization of the Hamiltonian is achieved, these two histogram should coincide with each 
other. Some theoretical argument on this issue is given in Appendix |A"1 

As we mentioned, this gap comes from the fact that the grayscale levels in the recognition image are restricted 
(defected) to s XyV — 0, 5, 10, 15, • • • , 255 due to the definition of the weight W x — i, y -j with size |N(x, y)\ <C 256. We 
discuss this issue later. In the right panel of FIG. |H we also show the histograms of grayscale levels for the original 
and threshold mask. Obviously, these two histograms are completely different as we also see the difference clearly 
from the upper left and the lower left panels in FIG. SJ 




From the view point of human visual systems, we check the power-spectrum of the threshold mask and the halftone 
binary dots. In order to investigate the statistical properties of the model, we next investigate the power spectrum of 
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the threshold mask and the halftone binary dots. Namely, we evaluate 

2 

\2n , „ . J 1 

: COS 



Po(fxJy) = { 72 £ ° l > 



L 2 ^ -J 



2tt 



L 2 ^ 



2vr 

-jr{ifx+3fy) 



(22) 



where we should replace Oij by the threshold mask tjj or the halftoned binary dots hij. We should bear in mind 
that here we consider the case L\ = L2 = L. 

We show the results in FIG. [7] As well-known, the so-called blue noise mask, there exists principal frequency 
and below the value, the power spectrum drops to zero, whereas, the high-frequency components remains finite. On 
the other hand, the threshold mask generated by our algorithm apparently depends on the structure of the original 
grayscale image. Therefore, it is assumed that some of the spatial structures of the original image, especially, the 
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FIG. 7: The power-spectrum of the halftone image (left) and threshold mask. 



smoothness leading up to the low frequency components in the power spectrum remains finite. However, from the 
power spectrum Pt{fx,f y ) shown in FIG[SJ we find that the low-frequency components almost disappear and high- 
frequency components remain finite. 

More convenient spectrum statistics is the so-called radial-averaged power spectrum 10] defined by 

Poifp) = N J (f)) E Po(fxJy), (23) 

where M(/ p ) means the region of annular rings in the Fourier space with the radiuses f p and /„ + A p , namely, the 
width of the annular rings is given by A p . iV(R(/ p )) is the number of frequency samples in i?(/ p ). 

In FIG. [HI we plot the spectrum statistics for both threshold mask and resulting binary dots of the halftoned image. 
Inset of FIG. [H shows the same statistics for the threshold mask image. From this figure, we find that the low-frequency 
components arc not dominant and relatively high-frequency components are larger than the low-frequency counter 
part. This tendency should be plausible from the view point of human-eyes modulation as we explained in FIG.[TJ 



IV. INVERSE DIGITAL-HALFTONING 



In the previous sections, we consider the digital halftoning from the view point of statistical mechanics of spin 
systems. In engineering perspective, it is also important for us to retrieve the original grayscale image from a given 
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FIG. 8: The radial-averaged power spectrum of the halftoned image (left) and the threshold mask (right). The inset in the left panel is 
the same plot as in the panel in the logarithmic scale. 

halftoned binary dots when we attempt to capture the halftone image via scanner machines. The inverse problem of 
the digital halftoning is referred to as inverse-halftoning. In this section, we investigate the statistical performance of 
the inverse-halftoning by making use of the concept of statistical mechanics of disordered spin systems. Especially, 
we discuss the condition on which the Bayes-optimal inverse-halftoning is achieved. 



A. Bayesian formula 



We first provide a Bayesian formula of inverse-halftoning. As we mentioned, the halftoned image h is given pixel- 
wise as h x>y — 0(g x ,y — t x ,y),^x,y Thus, we might regard the process as a kind of 'channel' which is written in terms 
of the following conditional probability: 



Pt(h\g) 



T\ x ,y S ( h x,y,®(9x,y ~ t x ,y)) 

Eh II , .., S(h x ,y, (-)'■!/, .., - t Xt y)) 



(24) 



for a given threshold mask t. Hence, the halftoning process is regarded as stationary memory-less deterministic erasure 
channel from the information theoretical point of view. When we choose the ferromagnetic prior, the posterior is 
obtained as 



P(tr\h : t) 



P t {h\a)Pj(a) 



exp[- J E((x, V )(fe,0) ^ X ' y ~ ak < 1 ^ 1LV !::/ ' ~ Q ^x,y - t x ,y)) 



E CT p t(h\cr)Pj(cr) E CT exp[- J £<(*,„)(*,,*)) (a x>y - cr w ) 2 ] U x , v 5 ( h x, y ~ ®{9x, y - t x , y )) 



(25) 



where we defined ^2/r xy \ (&;))(''') as a sum over the nearest neighbouring pixel pairs. From the view point of 
statistical mechanics, the above posterior is regarded as the Boltzmann distribution of the Q-Ising model in which 
the mobility of each spin is tightly constrained in the range t x y < g x y < Q — 1 for h x y = 1 and < g X/y < t x y for 



h X}V = 0. 



For the above posterior, the Maximizer of Posterior Marginal (MPM for short) estimate is given by pixel-wise 
optimization of the posterior as 



'x,y 



arg max P(o~ x , y | h : t) 



'E CT <fx,y e M-JY,{(x,y)(k,l)) ( a x,y ~ ^k,lf] Ux,y 5 ( h ^V ~ ®(9x,y ~ t x ,y)) ' 



E CT ex p[- J E 



(26) 



where Qq(x) stands for a function to convert a real value x to the nearest integer. 

In the inverse halftoning, we should solve V^.j, h x , v = Q(g x , y — t x .y) with respect to g for a given halftoned image 
h and the mask t. 
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In the reference [25|, one of the authors considered restoration processes of grayscale images by making use of 
bit-decomposed data. Namely, we generate the Q binary images whose pixel is given by V~ „ h^y = 0(gi — m), me 
{0, • • • , Q — 1} and transmitting them through some noise channels. Then, the paper |25| dealt with the procedure 
to retrieve the original image g from the degraded Q slices of the binary images 

h' (m) = h (m) +V m) (m = 0,l,--- ,Q- 1) (27) 

where we defined h = {hx$\x = 1, • • • , Li, y = 1, • • ■ ,L 2 } and r; (m) = {r)^y\x = 1, • • • , L\, y = 1, ■ ■ ■ ,L 2 } 
stands for the additive noise. If there is no degrading process, it is obvious that the vector 

O-i 

g = E hirn) ( 28 ) 

m=0 

is identical to the original grayscale image g. However, in the present inverse-halftoning case, the only information 
we have is just a single slice h^ m \ This fact makes the problem hard to treat. 

Obviously, this inverse-halftoning is a typical ill-posed problem because there are a lot of candidates to satisfy the 
equations. In Appendix [B] we evaluate two relevant quantities, namely, degree of degeneracy for possible solutions 
and mutual information to evaluate the difficulties of the problem. 

In the next subsection, we discuss the relationship between the so-called Bayes-optimal solution and the Nishimori 
line established in the research field of spin glasses. 



B. Bayes-optimal inverse-halftoning and the condition to achieve it 

In the previous studies [lH, we investigated the inverse-halftoning on the bases of Markov chain Monte Carlo 
simulations and analysis of the infinite-range mean-field model. Then, we found several conditions on the hyper- 
parameters appearing in the Hamiltonian that gives a minimum of the mean square error numerically. However, so 
far, we do not yet provide any mathematically rigorous results on the performance of the inverse-halftoning defined 
on realistic two dimensional square lattices. 

In this section, we attempt to prove that the Bayes-optimal inverse-halftoning is achieved on a specific condition 
which is similar to the so-called Nishimori line. For the purpose, let us use here an alternative definition of pixel index. 
Namely, for x,y = 1, 2, • • • , L, we change the index by means of i = x + Ly — L,i = 1, 2, • ■ • , L 2 = N (conversely, 
x = (i mod L), y — [i/L]). Then, we consider the true prior and the true inverse-halftoning process as follows. 

where we defined the original image g = {gi — 1, • • ■ ,Q\i — 1, • • ■ , N} and the halftoned image h — {hi = 0,l\i = 
0, • • • ,N}. We also used the definition (ij) to represent all nearest neighbouring pairs on the arbitrary lattice in 
finite dimension. The sums £ ff (- • • ) and E fc (---) denote £ g = n* E Si =i,- ,q0 ' ' ) and T, h = Hi Y, hi =o liG ' ' )> 
respectively. It should be noted that the above likelihood is another representation of the following dither method for 
each pixel 

V, hi = Qigi-ti). (30) 
For this original image and the halftone process given by the likelihood (|29p . we naturally use the following posterior: 

p , ... P t {h\<T)P jQ (<T) _ exp[-JE fa) (^-^) 2 ]n 8 ^ ; e(g«-*i)) 

The quantity to be evaluated is the following mean square error D for an arbitrary pixel i: 

d - E g , h[ { 9i - e Q ((a i}j ) } ] = E 9 ^M~JoE {v M-^mA^-u)) (32) 

where E g .h[- ■ •] = / J P t (h\g)Pj (g)(- ■ ■ )dgdh and we defined the average of a single pixel over the posterior by 
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To proceed the proof, we should notice that for any stochastic variable A, the following fact is satisfied E^jA — 
E^A]} 2 ] = E A L4 2 ] - E A [A} 2 > 0, namely, E A [A 2 } > E A [A} 2 . When we set A = g, - 6 Q ((cr 4 ). / ), we have the 
following inequality: 

D = E g , h [{g i -Q Q ((a i )j)} 2 ] > {E g , h [.9 J ]-E 9 ,h[e Q (( ( 7 4 ) J )]} 2 (34) 

In following, we shall evaluate the lower bound of the D, that is, the right most term of the above inequality. Among 
the staffs of the right most part in the above inequality, E 9 ,/i[<7i] is easily evaluated as 



E s,h[ft] = 



E g ,h ft exp[- Jp J2(ij) (ft - ft) 2 ] IL S(h z ,Q( gi - U)) 

J2 g ,h ex p[- J o E(ij) (ft - 9j) 2 \ Yli s{hi,e( gi - u)) 

Eg 9i exp[- Jp J2(ij) (ft ~ 9j) 2 ] Eh IX K h u Q(ft ~ U)) 
E 9 exp[- Jo E<ij) (ft ~ 9j) 2 ] Eh IL <K^> &(gi - U)) 

E. a «p[-As M ,( B -«,)■! _ (*r\ „,,,,><> + m „,,,,«> _ „,<« (36) 



where we defined 



Eg exp[- J E(y) (ft - ft) 2 ] V Zj 



m (Jo,9>t) = E ge{g ,| gi >t l} ft exph Jp Efa) (gi - ft) 2 ] 
E g e{ 9i | S4 >t 4 } ex Ph J o E<y) (ft - ft) 2 ] 
m (Jo, ff <t) = Egg {g ,| gi < tl} ft exp[- Jp ( gi - 9j ) 2 } 

Eg e{gi | 9l <t l} exp[- Jo £<y) (ft - 3,) 2 ] 



and 



4f 4) = E exp[-J ^( 5l - ffj ) 2 ] (38) 
ge{«|si>ti} <u) 

4o <4) = E exp[-JoE(ft-ft) 2 ] (39) 
ge{ g ,| gs <t s } (ii) 

Z Jo = ^exp[-JoE(ft-ft) 2 ] ( 4 °) 

9 <ii> 

Therefore, Eg.h[<7i] is identical to the local magnetization mf ^ of the pure ferromagnetic Q-Ising model having the 
interaction strength Jo. 

On the other hand, the lest of the term appearing in the right hand side of the equation (134)) is written by 
Eg,h[e Q ((a 4 ), 7 )] 

Eg,h eX P[" J 0E( li ) (ft - ft ) 2 ] Hi S {hh (ft - 

Eg exp[- Jo E(ij) (ft - ft) 2 ] 
/4f t} \ / E CTgKk ,> M ^exp[-JE fa) ^-^f] \ 
" 2j„ J ^ E CT{ . lkl >* l} exphJE fo) (a 4 -a,) 2 ] J 

( Z i <t] \ Qr (^e { a^<u } ^M-JE^(^-^) 2 }\ 



z Jo J V E CT { CT! | CTl<4l } ex Ph J E fo )(^-^) 2 ] 



• 7 (9>t)\ / 7 (g<t)\ 

^— \ e Q (m«">V ^- egf.f- 4 *) (4i) 



z 



Jo 7 \ ^Jo 
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Therefore, the lower bound of the mean square error is evaluated as follows. 
D > {E g , h [ 5l ]-E g , h [e Q ((a 4 ) J )]} 2 




where the equality on the last line is satisfied for J — Jo. We also used the fact that the sign of (mf°' 9> ' — 0Q(m/' <T<t )) 
and (mf°' 9> * — 0Q(wj/' ff< *)) are the same because the local magnetizations m J i ,a>t and mf a<t are monotonically 
increasing function with respect to J. 5 = S±, 62 is a quantization error due to the Q-generalized step function 0q(- • • ) 
defined as 

e Q (x) = x-S (43) 

From the argument we presented above, we found that the performance of the inverse-halftoning achieved by the 
posterior ()33p is optimized on the specific condition J — Jq which is similar to the Nishimori line (point) |2l| in the 
research field of spin glasses. 



V. SIMULTANEOUS GENERATION OF MASK, HALFTONE AND RECOGNITION IMAGES 

Finally we show that both halftoning and the inverse-halftoning processes are unified under a single Hamiltonian, 
namely, it is possible for us to obtain the threshold mask, the halftoned and inverse-halftoned images simultaneously 
by finding the ground state of the spin systems. 

In this paper, we proposed the digital halftoning and the inverse-halftoning separately. However, from the form 
of the Hamiltonian H(t\g), the definition of the halftoned image: V XlJ/ h XtV = &(g x , y — tx,y) and the recognition 
image: \/ XiV s x>y = J2i jeWx y) W x ~i,y-j®(gx,y — tx,y), it is confirmed that these three important images are obtained 
simultaneously in the single theoretical framework, that is, minimizing the Hamiltonian H(t\g) with respect to the 
mask t for a given original image g. 

Then, we should remember that there is a gap between the histograms of the grayscale levels for the recognition 
and the original images. This gap might make us hard to accept the recognition image as the inverse-halftoned one 
as a solution. 

As we mentioned before, this gap comes from the fact that the grayscale levels in the recognition image are restricted 
(defected) to s XlV = 0, 5, 10, 15, • • • ,255 due to the definition of the weight W x -i^ y ~j with size |N(x, y)| =5x5 = 
25 -C 256 = Q. To reduce the gap, we might use the following linear filter: 

k,l£N (x,y) 

where the initial condition are chosen as W x ,y s x . y — Yli jeM x y ) W x -i. y -jQ{g x . y — t x , y ) for the minimum energy 
state of the Hamiltonian 7i(t\g). Thus, for a given s^°> , we recursively operate the above map (|44p n times, and then, 
one might obtain more plausible image than the V X) « Sx°} y as the inverse-halftoned image. 

We show the resulting histogram of the grayscale levels in FIG. |9]and the corrected recognition images in FIG. [TO] 
for n = 1 and 2. We set the number of nearest neighbouring pixels around the pixel at (x, y) as |N (x, y)\ =3x3 = 9, 
namely, N'(x,y) = {(x-l,y + l), (x,y + l), (x+l,y+l), (x-l,y), (x,y), (x + l,y), (x-l,y-l), (x,y-l), (x+l,y-l)}. 
From these figures, we find that the gap between two images actually reduced and the resulting recognition image is 
improved by 'mixing effect' on the defected grayscale levels s x , y = 0, 5, 10, 15, • • • , 255. However, when we increase 
the number of iteration n, the mixing effect works too much on the recognition image and it makes the local structure 
of the image too smooth (see FIG. [S] (right) and FIG. [TU] (lower right)). 



y x,y o xy — 
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VI. SUMMARY 



In this paper, we proposed a possible statistical-mechanical modeling for digital halftoning. This formalism helps 
us to understand the problem as a combinatorial optimization which is described by a sort of disordered spin systems. 
Finding the ground state was achieved by simulated annealing and we found that the resulting binary dots looks 
well to represent the original grayscale levels. The quality of the binary dots was evaluated from the power-spectrum 
statistics. We found that the binary dots contain relatively high frequency components which are plausible from 
the view point of human-eyes modulation properties. We also proposed a theoretical framework to evaluate the 
statistical performance of the inverse digital-halftoning based on statistical mechanics. From the Bayesian inference 
view point, we rigorously show that the Bayes-optimal inverse-halftoning is achieved on a specific condition which is 
very similar to the so-called Nishimori line in the research field of spin glasses. We hope our formulation might be 
applied to generating of binary halftoned images effectively and evaluating the performance for the inverse-halftoning 
from halftoned binary images obtained by various algorithms. 

This work was financially supported by Grant-in-Aid, Scientific Research on Priority Areas Deepening and Ex- 
pansion of Statistical Mechanical Informatics (DEX-SMI) of the Ministry of Education, Culture, Sports, Science and 
Technology (MEXT) No. 18079001. One of the authors (JI) was financially supported by Grant-in-Aid for Scien- 
tific Research (C) of Japan Society for the Promotion of Science, No. 22500195 and INS A (Indian National Science 
Academy) - JSPS (Japan Society of Promotion of Science) Bilateral Exchange Programme. 



Appendix A: Distribution of the recognized pixels 



In Sec. IIIIl we evaluated the performance of halftoning through two different measurements, namely, the histograms 
of grayscale levels and the power-spectrum via computer simulations. However, it might be helpful for us to evaluate 
the performance analytically. In this Appendix, we derive the distribution of the recognized pixels, that is, the 
number of black pixels in a window with \N(x, y)\ — Q(-€i 256) analytically. Namely, we consider the distribution of 
the following quantities: 



1 



(Al) 



where we consider the case of L\ = 
i = 1, 2, ■ ■ • , I? and we set I? = N. 
we consider the case of — Q = 

ii ' 

Pp,Q(lH : flO = 



L2 = L and defined the index i so as to satisfy i = x + Ly — L, namely, i takes 
It is obvious that (a;, y) for a given i are given as x — (i mod L),y — [i/L]. Here 
: 0(1). Then, we obtain the general formula for the distribution as follows. 



-PH cl 



IL E^o I-oo IL diH8(QiH ~ Eiem)Q(gi ~ 



/3-Hcff 



(A2) 
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FIG. 10: From the upper left to the lower right, the original image, the recognized images after operating the linear filter d44[l with 
n = 1,2 and n = 20 are shown. 

where "Heff denotes the following effective Hamiltonian: 

H cS = (3'-^) + E E ®(9l-ti)®(9k-t k ) (A3) 

i ieN(i) i fe,/eN(i) 

In the above expression, we omitted the term 5? which is independent of the dynamical variable t. We should 
keep in mind that the normalization diiiP(fii : g) = 1 is satisfied. It also should be noted that the pixels are 
fully connected in the window with size |N(i)| = Q. Then, the effective Hamiltonian is reduced to the decoupled form 
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H e s = J2i ^cff an d we immediately obtain 



n e [°° n d ^ s f on e e ^ - *o ) ex p f - e ^ ) 

j tjl tt i J - 00 &i V !eN(i) / V j# / 

nE rn^M QMi ~ e (^-*o) exp(-^^F«) 

i ti i \ ZeN(i) J \ i J 



« - P -(W-1)/ 



M 1 - e -JV/ 



by using the saddle point method in the limit of iV — > oo. / is given by / = — $(/3, Q : g) with 

Q-l 



(A4) 
(A5) 



$(AQ:g) = log<j 
Then, equation (|A2|) leads to 



idp 
V2^ 



exp 



Qw* - f <2 V + + Q log ^ e 



fte(g-t) 



(A6) 



J — tea ^/2tt ^ 



-Qfi^ + Q log J2t=o 



exp 



Q-l 



-QAM - £Q V + + Q log ^ e 



AO(s-t) 



(A7) 



(A8) 



Obviously, the above distribution is dependent on the original image g and such data-averaged distribution is evaluated 
after slightly complicated algebra as follows. 



Pp, Q (fi) = Eg[P PiQ (f 1 :g)]=Eg 



h ex P (~ 2<3 2 ^ 2 + PQ9/j) f-oo dp cos [fip - Q tan 1 tp Q (p)] 



(A9) 



I- 00 d P z p.q(P) 

where Eg[- • • ] = J dg(- ■ ■ )P(g) and we defined the following functions. 

% Q (A) = &,q cos ^ iQ (A) + P (P) sin */3>q(A) + e"^ {^f ,Q) (£) cos - ,Q) (A) sin ^, Q (/i) JAIO) 



*/3,q(A) = Qtan" 



-i ( J2?Jo e^-t) smygAe(g - f) 
1 E?=o e^e( 9 -*) C os ^e(<? - t) 



^q(P) 

Po(P) 



E? = "o 1 cos(A/Q)e(. 9 -t) 



H(y/PQ), H{x) 



dt 



e 2 



(All) 

(A12) 
(A13) 



f 4^' ^ Q \P)^ r^e^cosVPQX, P^' Q HP) = f ^e^sin^QA (A14) 

JO V ^7T Jo V ^7T Jo V ^7T 



As a demonstration, we choose uniform grayscale images having a single grayscale go> namely, the distribution of G 
is given by P(g) = S(g — go)- Then, the data-average is easily performed as 



PpM ee E g [P((i:g)}=P(n:g ) 



(A15) 



In FIG. [Tl] we plot the P^ i q(/i) for several values of Q and /3 and g. From this figure, we find that as temperature 
fi^ 1 decreases, the distribution changes its shape to the delta function in which the peak is located at go/Q- This 
means that simulated annealing safely finds the ground state of the Hamiltonian having the zero energy, namely, 



Vj 5o = 9% = E ®( 5fc _ * fe ) = Si 



(A16) 



This immediately reads fa = go/Q for all index i. Therefore, for a given image having a single grayscale go, the 
distribution should converge to the delta function with the peak at go/Q if the annealing schedule of j3 is appropriate. 
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FIG. 11: The distribution P/j. q(/i) for several values of Q and /3 and g. 



Appendix B: Inverse halftoning as an ill-posed problem 



Here we shall consider the inverse process of halftoning as an ill-posed problem. In the inverse halftoning, we should 
solve y x ,y h x>y — 0(g x . y — tx,y) with respect to g for a given halftoned image h and the mask t. Obviously, this 
problem is ill-posed because there are a lot of candidates to satisfy the equations. In this Appendix, we analytically 
evaluate several relevant quantities to show the difficulty in finding the solution. 



1. Degree of degeneracy for possible solutions 

We easily assume that the number of the solutions for a given h and t is exponential order, however, it is helpful 
for us to evaluate the number more precisely. For the purpose, let us use the definition introduced in the previous 
section, namely, i = x + Ly — L,i = 1,2, •• ■ , L 2 = N (conversely, x = (i mod L), y = [i/L]). Then, the number of 
the candidates for the solution of the equations Vj hi — Q(gi — U) is evaluated for a given a realization of the original 
image and the threshold mask as follows. 

JV 

M(g>t) = n{(Q-l-ti)®(9i-ti)+U@(ti-gi)} (Bl) 

i=l 

As the number Af seems to be an exponential order object, we might rewrite the data average of the logarithm of the 
number as a self-averaging quantity 



JV 



log JJ {(Q - 1 - U)&{$i ~ U) + t&{ti - 9l )} 



(B2) 



E g 4\ogAf(g,t)} = E g ,t 
where we defined the expectation Eg,t[- • • ] by 

£„,*[■■■] = / [ P Jo (g)P(t)(---)dgdt, dg = Y[d gi ,dt = Y[dU. (B3) 



JV JV 

Y[d gi , dt=n 

" - i=l i=l 

Therefore, the above averages could be carried out for a specific choice of the P{g) and P(t). For distributions of the 
threshold mask, we suppose that in each mask with size Q x Q (Q <^ N), each component ti takes a value among 
0, • • ■ , Q — 1 grayscales with equal probability l/Q. Then, the P(t) is reduced to the product of the effective single 
site distribution as 

JV JV Q-l 

P(t) = ;Q p ^) = nllE^- ( B4 ) 

i=l ^ i=l U=0 
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On the other hand, as a distribution of the original image, we consider snapshots from the infinite-range ferromagnetic 
Q-Ising model, that is, 



Pjo(g) = 



exp(— H) 



2N 



(B5) 



where Zq is a normalization constant for the probability Pj (g). After simple algebra, we find that the Pj (g) is 



rewritten such as Pj (g) = Yli—i Pjoidi) with the following effective single site distribution: 

cxp(- J gf + 2J m Q g z ) T>a=o 9 exp(- J g 2 + 2J m g) 



Pj {m) 



to 



ElLo* ex P(- J o9 2 + 2J m g) 



(B6) 



where too = E-i=i 9i denotes the magnetization for the system of original grayscale images described by the 

Hamiltonian (|B5|) . For these probability distributions, we have the function <fi(Jo, Q) = A r_1 E g .t[logA/'(g, t)] explicitly 



as 



0(J o ,Q) = J2 Pj (9)P(t)log{(Q-l-t)Q(g-t) + tQ(t-g)} 



g,t=0 



E?=o E^o 1 kt E?=o ^ Jo9 2+2Jomo9 log {(Q - 1 - t)G(g - t) + tO(t g)} 



m g 



(B7) 



Therefore, the number of the solutions for the equations Vi hi — Q{g t — ti) is estimated as ~ exp(N 4>(Jq , Q)) . In 
FIG. H21 we plot the function (/>( J , 3) and 4>{ Jq, 4). From this figure, we find that the is a monotonically decreasing 



0.9 



4> 
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FIG. 12: Logarithm of the number of candidates (degeneracy) for the solutions <j> as a function of J 1 for Q = 3 and 4. 



function as J 1 decreases, however, the </> remains a finite value even if J 



d. This means that the number of the 



candidates for the solutions of the inverse halftoning is always exponential order. Therefore, we need some systematic 
approach to solve this type of the ill-posed problems. Namely, we need to introduce the ferromagnetic prior Pj(<t) to 
compensate the lack of information as we discussed in Sec. IIVI 



Mutual information 



We next consider the difficulties of retrieving original grayscale images from a slightly different point of view. Here 
we calculate the mutual information between the original image g and the halftone image h. From the definition 
of the mutual information, we should evaluate I(g : h) = H(h) — H(h\g), where we defined the entropy H(h) 
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and the conditional entropy H{h\g) as H(h) = -J2h P ( h ) ^ogP(h), H(h\g) = -£ fl PjM Eh p t(%) logP t (%), 
respectively. It should be noted that we used 

P t (h\g) = Hdih^&igi-ti)) (B8) 

i 

and P(h)=J2 p P t (h\g)Pj (g). 

For the infinite-range model and constant mask ti = to Vi, these entropies per are calculated analytically as 

H (h) ^ ^ f e- J °° 2 + 2J ° m ^8(h, Q(g - to)) \ , / E^o e- J ^ 2+2J ° m ^5{h, Q(g - to)) \ 

IT = -X>X\ zj los ^ < ( } 



H(h\g) ^ e 



h=0,l 9=0 

1 ( -J g 2 +2jQm g 



v E — z r 5 ^ e (9~ <o )) lo § e ^ - *<>)) = ( B1 °) 



i — ' Jn i 

9=0 k Jo j h=01 



where we defined Zj = E^Lo* e J " 9 +2Jom " 9 and mo is a solution of the equation (|B6I) . We should keep in mind that 
the 'channel' of the halftoning process is completely deterministic. As the result, the conditional entropy H(h\g) is 
identically zero. In FIG. [131 we plot the mutual information per pixel I(g : h)/N as a function of Jq. To compare 



Halftone 
aja=0.5 
aja=l 
a ! Ja=2 




FIG. 13: The mutual information per pixel as a function of Jo- For the halftoning process, we choose Q = 3 and to = 1. The mutual 
information for image restoration under Gaussian noise is also plotted for the signal to noise ratio an. /a = 0.5, 1 and 2. 

halftoning process by the dither method with the degraded process by a Gaussian noise with mean a^g and the variance 
a 2 , we also calculate the mutual information per pixel for the Gaussian noise. We choose the same distribution Pj (g) 
as the halftone case. We immediately have 

H(h) J f dh y e - Jog * +2Jomog -2^1- lQg U=o e (Bn) 



N V2iraZj J_oo ^ " I V2iraZj 

Wi± = i+log(V<M. (B12) 

In FIG. [131 we plot the mutual information per pixel for the Gaussian noise as a function of Jo for various cases of 
the signal- noise ratio, Qq /a. 
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